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ABSTRACT 

We investigate the gravitational interaction of a Jovian mass protoplanet with a 
gaseous disc with aspect ratio and kinematic viscosity expected for the protoplan- 
etary disc from which it formed. Different disc surface density distributions have been 
investigated. We focus on the tidal interaction with the disc with the consequent gap 
formation and orbital migration of the protoplanet. Nonlinear hydrodynamic simula- 
tions are employed using three independent numerical codes. 

A principal result is that the direction of the orbital migration is always inwards 
and such that the protoplanet reaches the central star in a near circular orbit after a 
characteristic viscous time scale of ~ 10 4 initial orbital periods. This was found to be 
independent of whether the protoplanet was allowed to accrete mass or not. Inward 
migration is helped through the disappearance of the inner disc, and therefore the 
positive torque it would exert, because of accretion onto the central star. Maximally 
accreting protoplanets reached about four Jovian masses on reaching the neighbour- 
hood of the central star. Our results indicate that a realistic upper limit for the masses 
of closely orbiting giant planets is ~ 5 Jupiter masses, because of the reduced accretion 
rates obtained for planets of increasing mass. 

Assuming some process such as termination of the inner disc through a magne- 
tospheric cavity stops the migration, the range of masses estimated for a number of 
close orbiting giant planets (Marcy, Cochran, & Mayor 1999; Marcy & Butler 1998) 
as well as their inward orbital migration can be accounted for by consideration of 
disc-protoplanet interactions during the late stages of giant planet formation. 

Key words: giant planet formation- extrasolar planets- accretion discs- numerical 
simulations 



1 INTRODUCTION 

The recent discovery of a number of extrasolar giant plan- 
ets orbiting around nearby solar-type stars has stimulated 
renewed interest in the theory of planet formation. These 
planetary objects have masses, m p , that are comparable to 
that of Jupiter (0.4 Mj < m p < 11 Mj), have orbital semi- 
major axes in the range 0.04 AU < a < 2.5 AU, and orbital 
eccentricities in the range 0.0 < e < 0.67 (Marcy, Cochran, 
& Mayor 1999, Marcy & Butler 1998 and references therein). 
It should be noted that the detection technique of measur- 
ing the Doppler shift induced by the host star's orbital reflex 
motion only allows the measurement of m p sini, where i is 
the inclination angle of the orbit plane to the line of sight. 
It is generally believed that planets form out of the gas 
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and dust contained in the discs that are observed around 
young T Tauri stars (Beckwith & Sargent 1996). This T 
Tauri disc phase of a stars life is thought to last on the order 
of 10 6 -10 7 yr, after which the discs appear to dissipate. In 
the standard theory, planet formation occurs through a num- 
ber of key stages. First, the dust grains, that are initially well 
mixed with the gas in the disc, undergo coagulative growth 
via binary collisions. Second, as the grains continue to grow 
they begin to gravitationally settle towards the midplane 
of the protostellar disc, forming a dense dust layer in the 
process. The existence of this dense layer enhances the rate 
at which the solid material may combine into larger bodies, 
leading eventually to the formation of planetesimals. Third, 
the planetesimals continue to grow through collisions, pos- 
sibly aided by a runaway accretion process (e.g. Lissauer & 
Stewart 1993) ultimately forming planetary sized objects. 
These authors estimate that the time scale for this to occur 
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at 5 AU is 10 5 -10 6 yr, although this would require a high 
dust to gas ratio in a minimum mass solar nebula. However, 
there are many uncertainties in the processes involved and 
the effect of phenomena such as disc-planet interactions and 
orbital migration considered in this paper have yet to be ex- 
plored in detail. The latter might decrease the time scale by 
providing increased mobility of protoplanets in the nebula. 

In the 'critical core mass' model of giant planet for- 
mation, the formation proceeds through the build-up of a 
critical ~ 15 M eart h solid core, beyond which mass rapid 
gas accretion occurs leading to the formation of a gas giant 
planet. Because this process must occur prior to the loss of 
the gas from the disc (i.e. within ~ 10 7 yr) it is expected 
that the cores of gas giant planets should begin to form be- 
yond a radius of r ~ 4 AU, the so-called 'ice condensation 
radius' or 'snow-line'. Beyond this radius the temperature in 
the protostellar disc falls below the level that allows volatile 
materials to condense out into the solid phase and form ices. 
The presence of this additional solid material increases the 
rate at which solid materials may coagulate into larger ob- 
jects, and thus shortens the time required to form the solid 
cores that are thought to be the precursors to gas giant 
planets. 

If giant planets begin to form beyond radii of ~ 4 AU 
from their host stars, then orbital migration must have oc- 
curred in order to explain the existence of the closely or- 
biting extrasolar giant planets. The most likely cause of 
this orbital migration is through the gravitational interac- 
tion between the embedded protoplanet and the protostellar 
disc. The linearised response of a gaseous disc to the pres- 
ence of an embedded satellite has been investigated by a 
number of authors (e.g. Goldreich & Tremaine 1978, 1979, 
1980; Lin & Papaloizou 1979a,b, 1980, 1993; Papaloizou & 
Lin 1984; Ward 1997 and references therein; Artymowicz 
1993a, b). The perturbation of an accretion disc by a proto- 
planet leads to the excitation of spiral density waves at Lind- 
blad resonances which carry with them an associated angu- 
lar momentum flux that is deposited in the disc at the loca- 
tion where the waves are damped. The disc orbiting beyond 
the position of the planet receives angular momentum from 
the planet, whereas the inner disc loses angular momentum 
to the planet. In the situation where the tidal torques are 
greater than the internal viscous torques in the disc, and 
the disc response becomes non linear, it is expected that an 
annular gap, or surface density depression, may be formed 
in the vicinity of the planet (Papaloizou & Lin 1984, Lin & 
Papaloizou 1993). This tidal truncation of the protostellar 
nebula was investigated using non linear numerical simula- 
tions by Bryden et al. (1999), hereafter BCLNP, and Kley 
(1999a), in order to examine the effect of gap formation on 
the mass accretion rate by an embedded giant protoplanet. 
The results of these studies indicate that, for physical pa- 
rameters typical of protostellar disc models, gap formation 
can substantially reduce the accretion rate, leading to ex- 
pected planet masses in the range 1 < m p < 10 Mj, in 
close agreement with the observed masses of the extrasolar 
planets. 

The exchange of angular momentum between the planet 
and the disc leads to the possibility of the planet undergo- 
ing orbital migration if there exists an imbalance between 
the torques exerted by the inner and outer discs (Goldre- 
ich & Tremaine 1980). The migration that results may be 



described by one of two different formalisms, depending on 
whether the disc response is linear or nonlinear (with a gap 
forming) . 

Type I migration occurs when the disc response is lin- 
ear and the background surface density profile remains es- 
sentially unchanged by the interaction. A natural tendency 
for the outer disc torques to dominate the inner disc torques 
results in the inward migration of the planet independently 
of the background disc flow, and typically results in rapid 
migration time scales of 10 5 or 10 4 yr for a 1 or 10 M ear th 
planet located at 5 AU, respectively (Ward 1997). 

Type II migration occurs when the disc response be- 
comes non linear and a clear gap is formed around the vicin- 
ity of the planet. In this case, provided the planet mass is 
less than or comparable to the local disc mass with which 
it interacts, the migration occurs on the viscous evolution 
time scale of the disc. This process was investigated in detail 
by Lin & Papaloizou (1986). When the mass of the planet 
becomes larger than the local disc mass, then the inertia 
of the planet becomes important in determining the migra- 
tion rate of the planet. This situation has been investigated 
by Syer & Clarke (1995), and more recently by Ivanov, Pa- 
paloizou & Polnarev (1999) in the context of satellite black 
holes orbiting in AGN discs. 

Both type I and type II migration may be important at 
different stages of the planet formation process. In particu- 
lar it is possible that solid cores undergo type I migration 
before gas accretion leading to giant planet formation at 
radii smaller than 4 AU (Papaloizou & Terquem 1999) and 
possibly in situ close to the central star. 

Extrasolar giant planets, however, are observed to orbit 
over a wide range of radii, some of which must have been 
contained within the original gas disc prior to its dissipation. 
Accordingly, in this paper we shall focus on the situation 
where a giant planet is assumed to have formed at some 
radius in the disc and investigate the subsequent type II 
migration that occurs as a consequence of disc-protoplanet 
interaction. 

We employ three independent Eulerian hydrodynamic 
codes to examine the non linear evolution of a combined star, 
disc, and embedded protoplanet system. The main questions 
that we wish to address are: 

(i) what is the time scale for an initially embedded giant 
protoplanet to migrate in to the close vicinity of its host 
star ? 

(ii) How much mass is the planet able to accrete from the 
disc during this time, and how does the mass evolution of 
the planet affect its orbital migration rate ? 

(iii) Does the tidal interaction between the disc and 
planet cause the growth of orbital eccentricity during the 
migration ? 

The physical parameters that we focus on in this paper 
are appropriate to a 1 Mj protoplanet embedded in a min- 
imum mass solar nebula model containing ~ 2 Mj within 
5 AU (i.e. the initial orbital radius of the protoplanet). We 
find that the protoplanet migrates in towards the host star 
approximately on the viscous evolution time of the disc, in- 
dependently of the details of the initial conditions of the 
simulations, or the numerical code used. For a planet ini- 
tially located at a radius of 5 AU this translates into a time 
of ~ 10 5 yr for the disc parameters that we employ. This 
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time is substantially shorter than the expected life times of 
protoplanetary discs, and indicates that orbital migration is 
an important factor during the formation epoch of all plan- 
etary systems. 

Assuming maximal accretion, the estimated final 
masses of the planets as they approach their host stars is 
found to be in the range 2.7 < m p < 5.65 Mj, depending on 
the details of the calculation. These values fit in well with 
the observed mass range of the extrasolar planets. 

The orbits of the planets in all calculations were found 
to remain essentially circular. Thus the observed eccentric- 
ities of the extrasolar planets are not reproduced by our 
current models. These might be produced if the disc had 
a lower viscosity resulting in wider and deeper gaps than 
obtained here or by the perturbing presence of additional 
planets in the system. These issues will be the subject of a 
future investigation. 

This paper is organised as follows. In section || we 
present a more quantitative discussion of gap formation and 
orbital migration. In section ^ we present a brief discussion 
of protostellar disc models, and our choice of disc param- 
eters. This is followed by a discussion of the equations of 
motion, boundary and initial conditions, physical parame- 
ters, and treatment of the protoplanet in section ^[ We then 
go on to describe the hydrodynamic codes that we use in 
. The results of the calculations are described in 
There, we discuss in detail the results of one cal- 
culation (our standard run), and then examine how the re- 
sults depend on: (1). the presence or not of an initial gap in 
the vicinity of the planet; (2). whether the planet is accret- 
ing or non accreting; (3). the numerical resolution; (4). the 
code used to perform the calculation. We also present the 
results of one very long time scale evolution calculation and 
investigate the effects of changing the initial density profile. 
Finally, in section (7) we discuss the broader implications of 
our results and draw our conclusions. 



2 ORBITAL MIGRATION AND GAP 
FORMATION 

The tidal interaction between an accretion disc and an em- 
bedded protoplanet leads to the exchange of angular mo- 
mentum between them. The nonaxisymmetric surface den- 
sity perturbation of the more slowly rotating disc exterior 
to the orbital radius of the planet produces a negative grav- 
itational torque acting on the planet. Similarly the more 
rapidly rotating inner disc exerts a positive torque. Any im- 
balance between these inner and outer torques will lead to 
the orbital migration of the planet (Goldreich & Tremaine 
1980). 

By Newton's third law the planet exerts oppositely di- 
rected torques on the inner and outer disc material. An an- 
nular gap may be formed locally in the disc if the magni- 
tudes of these torques exceeds the internal viscous torques 
(Papaloizou & Lin 1984). The presence or absence of a gap 
determines whether the migration is of type II or I, respec- 
tively (Ward 1997). 



2.1 Type I Migration 

Type I migration occurs when the response of the disc to 
tidal forcing by an embedded planet is linear. Then no gap 
is formed and the background surface density profile is ap- 
proximately unchanged. 

The presence of a planet orbiting in a gaseous disc leads 
to the excitation of trailing spiral density waves at the Lind- 
blad resonances in the disc (Goldreich & Tremaine 1979). 
These density waves carry with them an associated angular 
momentum flux which is deposited locally in the disc at the 
location where the waves damp. The gravitational coupling 
between the trailing wave pattern and the planet leads to a 
torque acting on the planet. 

A natural imbalance between the torques acting on the 
outer and inner disc arises because the locations of the outer 
Lindblad resonances tend to be closer to the planet's posi- 
tion than the inner ones, leading to a net inward migration 
of the planet. This sense of migration is insensitive to details 
of the background disc flow (e.g. Ward 1997). 

The differential torque induced by the Lindblad reso- 
nances may be written as (Ward 1997) 



AT 



(|)^(^) 2 (^) ; 



(1) 



where q is the planet to central star mass ratio m p /M», E 
is the surface density of the disc material, r v is the planet 
orbital radius, £1 is the Keplerian angular velocity, H/r is the 
disc aspect ratio, and ci accounts for the torque imbalance 
between the two sides of the disc and is expected to scale as 
~ (H/r). The corresponding radial velocity of migration of 
the planet is given by 



dr p 
dt 



-aq 



M, 



(2) 



This equation predicts that the radial migration will occur at 
a higher rate for a more massive planet, and remains valid 
until the disc response becomes non linear and a gap be- 
gins to form. A sufficient condition for non linearity through 
shock formation is (e.g. Korycansky & Papaloizou 1996) 



q > 



(3) 



This corresponds to a planet of mass ~ 30 M eart |j in a proto- 
stellar disc with a typically expected aspect ratio H/r ~ 0.05 
(see Papaloizou & Terquem 1999 and references therein). 
Accordingly we expect, as is confirmed by the simulations 
presented here, that type I migration applies to protoplanets 
of substantially smaller mass than we consider in this paper. 



2.2 Type II Migration 

When the disc response becomes non linear a gap forms, 
inside of which the planet orbits. The conditions for a gap 
to form have been discussed by Papaloizou & Lin (1984), 
Lin & Papaloizou (1985, 1993) and BCLNP. These are the 
thermal or shock formation condition given by Eq. ([]) and 
the condition that tidal torques exceed viscous torques which 
may be written 
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where 1Z = r 2 Q,/v is the Reynolds' number, and v is the 
kinematic viscosity. 

When there is a gap and the planet mass is less than or 
comparable to the local disc mass with which it gravitation- 
ally interacts, the migration rate of the planet is controlled 
by the viscous evolution of the disc, since the planet then 
behaves as a representative particle in the disc. In this case 
the migration rate is given by the radial drift velocity of the 
gas due to viscous evolution, which for a steady state disc is 
given by 

dr v ?>v 

(5) 



dt 



2?' D 



This leads to a migration time 
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When the mass of the planet is larger than the charac- 
teristic disc mass with which it tidally interacts, the inertia 
of the planet becomes important in slowing down the orbital 
evolution. The inertia of the planet acts as a dam against 
the viscous evolution of the disc, and can lead to a sub- 
stantial change in the disc structure in the vicinity of the 
planet. The coupled disc-planet evolution in this case has 
been considered by Syer & Clarke (1995) and more recently 
by Ivanov, Papaloizou & Polnarev (1999). Using a simple an- 
alytical model, Ivanov el al. (1999) estimate the migration 
time of a massive planet to be 

1/5 



1 M. 



,/n 



10 m„ 



M 



for a disc with constant v, where M is the mass accretion rate 
through the disc and Mdo is the characteristic unperturbed 
disc mass that would be contained within the orbital radius. 
We can write M = Mdo/i~ v (r p ), where r v (r p ) — 2r p /(3v) is 
the viscous evolution time of the disc at a distance r p , and 
Mdo = leading to 

1/5 
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If we write r m i S ~ r p (dr p /dt) 1 , then we obtain the follow- 
ing relation for the migration rate 



^ ~ f (10^E 4 ) 1/5 m p - 4/5 rr. 



(9) 



Thus, we see that as the mass of the planet increases, or its 
orbital radius decreases, the rate of migration should slow 
down. The latter effect arises because the planet interacts 
with a smaller amount of disc mass at smaller radii. This 
analysis also predicts that a protoplanet with mass substan- 
tially larger than Mdo should not increase its mass signifi- 
cantly before migrating to the centre of the disc. This to- 
gether with the reduction of the accretion rate with increas- 
ing protoplanet mass (BCLNP and work presented here) 
suggests that the protoplanet mass should also be limited 
at about Mdo- 

If we consider the interaction of a Jupiter mass planet 
initially at 5 AU with a minimum mass solar nebula disc 
model containing ~ 2 Jupiter masses within 5 AU, we see 
that the migration rate should initially occur at the viscous 
evolution rate of the disc given by Eq. Q since m p < Mdo- 
However, if the planet accretes mass and/or migrates in- 
wards, then m p eventually becomes larger than Mdo, such 



that Eq. (|9j) may apply. Thus, the parameter regime that we 
consider in this paper is expected to be transitional between 
those governed by Eqs. (0) and (H). 



3 PROTOSTELLAR DISC MODELS 

Models of protostellar discs considered as viscous accretion 
discs have been constructed by a number of authors (e.g. Bell 
el al. 1997, Papaloizou, Terquem, & Nelson 1999). An impor- 
tant issue is the nature of the effective viscosity. Usually the 
'a' prescription of Shakura & Sunyaev (1973), v — aH 2 Q is 
adopted. The most likely source of the turbulence required to 
produce the effective viscosity is MHD instabilities (Balbus 
& Hawley 1991, 1998). Simulations have shown that values 
of a in the range 10" 3 -1(T 2 may be produced. However, it 
is unclear that adequate ionization exists for the mechanism 
to work throughout the disc and in some cases MHD insta- 
bilities may exist only in a surface layer ionized by cosmic 
rays (Gammie 1996). 

Assuming MHD instabilities do work throughout the 
disc and produce values of a in the above range, disc models 
with properties similar to that of a minimum mass solar 
nebula are produced at a time ~ 10 5 - 10 6 yr after formation 
(Papaloizou, Terquem, & Nelson 1999). These typically have 
H/r ~ 0.04 — 0.05 so that the kinematic viscosity v ~ 10~ 5 
at 5 AU. 

Although there are uncertainties as to how a solid core 
of sufficient mass accumulates for rapid gas accretion to be- 
gin, in order to understand the distribution of the orbital 
parameters of extrasolar planets, it is reasonable to pose the 
question as to how a Jupiter mass protoplanet evolves as a 
result of interaction with the protostellar disc immediately 
post formation. 

Accordingly we shall consider the interaction of a 
Jupiter mass protoplanet with a disc containing two Jupiter 
masses interior to the initial protoplanet orbit and with a 
constant v = 10~ 5 . We consider cases where the initial disc 
surface density is uniform and where it depends on radius. 



4 THE PHYSICAL MODEL 
4.1 Equations of Motion 

The vertical thickness H of an accretion disc, which is in a 
state of near-Keplerian rotation, is small in comparison to 
the distance r from the centre, i.e. H/r <C 1. It is therefore 
convenient to vertically average the equations of motion and 
work with vertically averaged quantities only, the assump- 
tion being that there is zero vertical motion. The problem 
is thus reduced to a 2-dimensional one. 

We shall work with cylindrical coordinates (r, <f), z), 
with the origin located at the position of the central star. 
The velocity in our 2-dimensional disc is denoted by v = (v r , 
1)0, 0), where v r is the radial velocity and is the azimuthal 
velocity. We denote the angular velocity of the disc material 
by S7 = w^/r, where the rotation axis is assumed to be coin- 
cident with the vertical axis of the coordinate system. The 
vertically integrated continuity equation is given by 

9E 
dt 



+ V.(Sv) =0. 



(10) 
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The components of the momentum equation are 

dCSVr) . Vol gp , 

_L^> + v . (EwrV) = _^ _ _ E + fr 



(11) 



9t 



+ V.(£tuv) 



-^-^--1? + /, (12) 



Here E denotes the surface density 

/oo 
pdz, 
- oo 

with p being the density, P is the vertically integrated pres- 
sure, and f r and f# denote the viscous force per unit area 
acting in the r and 6 directions respectively. The gravita- 
tional potential, $, is given by 

GM, Gm„ 



y/r 2 + rl 



p) 



+ 



Gm v 



G 



dm(r') 



(13) 



where M* and m p are the masses of the central star and 
the protoplanet respectively, and r v and 6 V are l ne radial 
and azimuthal coordinates of the protoplanet. The third and 
fourth terms in Eq. ( |l3|) account for the fact that the coor- 
dinate system based on the central star is accelerated by 
the combined effects of the protoplanet ary companion and 
by the gravitational force due to the disc, respectively, with 
the integral in Eq. ([[3]) being performed over the surface of 
the disc. 

In our models the protoplanet evolves under the grav- 
itational attraction of the central star and the protostellar 
disc. The latter interaction is expected to cause the proto- 
planet to undergo orbital evolution. The equation of motion 
of the protoplanet may be written 



d 2 r p _ G(M« + mp) 
~~dlF ~ rl 



(14) 



where the gravitational potential of the disc is given by 



$ d = -G 



|r' - r„ 



■dr' + G 



dm(r') 



(15) 



Here the integrals are performed over the disc surface, and 
the second term is the indirect term arising from the fact 
that the coordinate system is accelerated by the disc grav- 
ity — note that the part of the indirect term due to the 
planet itself is already included in Eq ([14]). As the proto- 
planet orbits about the central star, it is able to accrete gas 
from the surrounding protostellar disc. The disc gas that 
it accretes has an associated specific angular momentum, 
which if different from the specific angular momentum of 
the protoplanet, will cause its orbit to evolve. We include 
the effects of this 'advected angular momentum' on the pro- 
toplanet's orbit by calculating the appropriate changes to 
the planet's mass and angular momentum that arise from 
the gas accretion at each time step. 



4.2 Equation of State 

For computational simplicity we adopt a locally isothermal 
equation of state. The vertically integrated pressure is re- 
lated to the surface density through the expression 



where the local isothermal sound speed is given by 

H 

Cs = VK 

r 



(16) 



where vk = \J GM, / r denotes the Keplerian velocity of the 
unperturbed disc. The disc aspect ratio, H/r, is assumed 
to be an input parameter that defines the Mach number of 
the disc model being considered. The calculations presented 
in this paper for the most part adopted H/r — 0.04. Some 
calculations denoted by Ri adopted H/r = 0.05, see table |l| 

4.3 Viscosity 

In this present work, we assume that protostellar discs have 
an anomalous effective viscosity most probably arising from 
MHD turbulence (see discussion in section ^). Here we make 
the assumption that this effective viscosity can be modeled 
by simply replacing the molecular kinematic viscosity coeffi- 
cient in the Navier-Stokes equations by a turbulent viscosity 
coefficient denoted by v. 

The components of the viscous force per unit area may 
then be written 



_ 1 d(rr r , 
r dr 

1 <9(rT0 r ) 

U = a 

r or 



1 dr r ^ 
r 86 



T <t>4> 
r 



1 dr^ Trj, 
r dd> r 



(17) 



(18) 



where the components of the viscous stress tensor used in 
the above expressions are 



2rjD rr 



2 „ 



(19) 



2r)D^ - -77V.V 



where 



Drd> — r 



dv r 
dr 
1 
2 



D<f>4> = — 



1 dv<t, 



+ — 

r 



(20) 



.JL (^±) + l dVr 

dr \ r J r 86 



and rj = Ez/ is the vertically integrated dynamical viscosity 
coefficient. In the work presented in this paper we use a 
constant value for v = 10 -5 . 



4.4 Dimensionless Units 

For reasons of computational simplicity, we use dimension- 
less units for our numerical calculations. The unit of mass is 
taken to be the sum of mass of the central star M* and the 
initial mass of the protoplanet m p , and the unit of length is 
taken to be the initial orbital radius of the protoplanet 00. 
We set the gravitational constant G = 1. The unit of time 
then becomes 



to 



G(M* + m p ) 



When discussing the results of the calculations in subsequent 
sections, we will use the initial orbital period of the planet 
as the unit of time, given by Pq = 2-nto. 
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4.5 Gas Accretion by the Protoplanet 

Following Kley (1999a), the accretion of gas by the proto- 
planet is modeled by removing a fraction of the material 
that resides within a distance of |i? r from the protoplanet 
at each time step, and a different fraction from within \R r , 
where R r is the Roche radius of the planet given by 

^G^r (21) 
The fraction of gas removed at each time step determines the 
local accretion time scale onto the protoplanet, t acc , which 
in our calculations is taken to be t aC c = 3to within l-R r and 
tacc = l-5to within jR r - This accretion rate is large and 
is almost maximal (Kley 1999a). We also perform simula- 
tions in which the accretion rate is set to zero. In this case 
a lobe filling atmosphere develops such that material that 
approaches the protoplanet is forced by its pressure to either 
return or cross the gap. Thus our calculations span the range 
of possible accretion rate behaviour onto the protoplanet. 

4.6 Initial Conditions 

The disc models used for all simulations using the codes 
NIRVANA and FARGO were of uniform surface density ini- 
tially, had a value of v = 10 -5 throughout and a constant 
value of H/r = 0.04. The value of E was chosen such that 
there exists the equivalent of 2 Jupiter masses in the disc 
interior to the orbital radius of the protoplanet, initially. In 
our dimensionless units this gives E = 6 x 10~ 4 . The ini- 
tial mass ratio between the protoplanet and the central star 
was taken to be q = m p /M* = 10 -3 , corresponding to a 
Jupiter mass planet orbiting about a solar mass star. The 
protoplanet was started on a circular orbit of radius r = 1. 
The inner radius of the disc is located at r = 0.4 and the 
outer radius is located at r = 2.5. 

Simulations with RH2D were carried out using identical 
initial conditions to those described in Kley (1999a). The 
aspect ratio of the disc was H/r = 0.05. The initial density 
profile was of the form E(r) oc r~ 1//2 , and in all cases an 
initial annular gap was imposed. The inner disc radius is 
located at r = 0.25 and the outer radius is at r = 4.0. . 
The total initial disc mass was 10 _2 Mq. In comparison, the 
initial models for NIRVANA and FARGO have a mass 1.45 
times larger. Hence, the surface density of the models with 
RH2D at r = 1 is about three times lower. 

4.7 Boundary Conditions 

An outflow boundary condition is used at the inner bound- 
ary for all calculations presented here since material in a 
viscous accretion disc will naturally flow onto the central 
star. 

The outer boundary condition is more problematic, 
since ideally we would like to have a closed outer boundary. 
We work in a coordinate system that is based on the central 
star and not on the centre of mass. The natural tendency for 
material at the outer edge of the disc is to orbit about the 
centre of mass and not the central star. Adopti ng the u sual 
closed boundary condition of v r = and = y/ GM„ / r will 
result in a mismatch between the disc material just interior 
to the outer boundary and that imposed at the boundary 



itself. Calculations using this condition do indeed show the 
resulting excitation of waves at the outer edge of the disc, 
even for the case of a non-accreting planet on a fixed circular 
orbit. 

In order to alleviate this problem, the boundary condi- 
tion that we adopt in the NIRVANA calculations assumes 
that material at the outer edge of the disc is in circular, 
Keplerian orbit about the centre of mass of the star plus 
planet system. Since we work in a coordinate system based 
on the central star, this requires us to calculate the correct 
values of v r and «0 at each time step and apply them to the 
outer boundary. One consequence of this is that the outer 
boundary is no longer closed, but can allow the inflow and 
outflow of material since the radial velocity is no longer zero. 
This effect is small, however, and appears to have a negli- 
gible effect on our simulation results. Simulations of a non- 
accreting planet on a fixed circular orbit using this boundary 
condition show no signs of wave excitation at the outer edge 
of the disc. When the planet is able to accrete gas from the 
disc and to orbitally migrate, this boundary condition again 
results in some wave excitation at the outer edge, but with 
a much reduced amplitude. This arises because the centre 
of mass position changes with time as the planet grows in 
mass and changes its orbit, and these changes are fed into 
the outer boundary condition instantaneously. The disc ma- 
terial interior to the outer boundary on the other hand, will 
adjust to the centre of mass evolution on a longer time scale, 
so that this boundary condition also produces a small but 
noticeable mismatch between the boundary and the outer 
disc material. 

The calculations performed using FARGO and RH2D 
all used a closed outer boundary condition such that v r = 0. 
The strong similarities in the results of the calculations per- 
formed with NIRVANA and FARGO indicate that the de- 
tails of the outer boundary have a negligible effect on our 
results. This is because, as tests have shown, the protoplanet 
primarily interacts with material that is close to its imme- 
diate vicinity. 



5 THE HYDRODYNAMIC CODES 

In order to establish the reliability of the numerical results, 
the equations of motion ( |io| ) to (|l^) have been solved using 
three different Eulerian hydrodynamic codes, NIRVANA, 
FARGO and RH2D. In each case the equations are solved 
using a finite difference scheme on a discretised computa- 
tional domain containing N r x N$ grid cells. Each scheme 
is described briefly below. 

5.1 NIRVANA 

NIRVANA is a three dimensional MHD code that has been 
described in depth elsewhere (Ziegler & Yorke 1997). For the 
simulations presented here, the magnetic field is set to zero 
such that the code becomes purely hydrodynamic. We work 
in two dimensions and use cylindrical (r, <j>) coordinates. 
Viscous terms have been added as described by Kley (1998). 
The computational domain is subdivided into NrXN^ zones, 
where the grid spacing in both coordinate directions is uni- 
form. 

For the calculations that are presented in this paper, 
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Figure 1. This figure shows the evolution of a protoplanct embedded in a protostellar disc for calculation N3. The relative surface 
density of disc material is represented by the grey-scale. The white circle represents the position of the protoplanct. The disc is initially 
unperturbed at time t = 0. 



three different levels of resolution have been used. The low 
resolution runs use iV r = 50 and JV^ = 147, the mid- 
resolution runs use iV r = 80 and — 235, and the high 
resolution runs use N r — 130 and N^, — 384. The numer- 
ical method is based on a spatially second-order accurate, 
explicit method that computes the advection using the sec- 
ond order monotonic transport algorithm (Van Leer 1977), 



leading to the global conservation of mass and angular mo- 
mentum. The evolution of the planet orbit is computed using 
a standard leapfrog integrator. NIRVANA has been applied 
to a number of different problems including that of an ac- 
creting protoplanet embedded in a protostellar disc. It was 
found to give results that are very similar to those obtained 
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with other finite difference codes including RH2D (e.g. Kley 
1999a). 

5.2 FARGO 

This is an alternative Eulerian ZEUS-like code, based on the 
FARGO fast advection method (Masset 1999). The main dif- 
ference between NIRVANA and this code is that in FARGO 
the time step is not limited by the classical CFL condition, 
which results in a very small time step due to the fast orbital 
motion at the inner boundary, but it is limited by a CFL 
condition based on the residual velocity with respect to the 
average orbital motion. This leads to a substantially larger 
time step and hence faster computation. Also, the FARGO 
procedure leads to a smaller numerical diffusivity because a 
larger time step size requires one to perform fewer advec- 
tion substeps during the calculations. Since the time step in 
the FARGO simulations can be quite large, especially in the 
low resolution case, a fourth-order Runge Kutta scheme was 
used to integrate the equations of motion of the protoplanet. 

Since this work represents the first application of the 
FARGO advection algorithm, it has been widely tested 
against NIRVANA by running strictly similar simulations 
(identical physical and numerical parameters). The good 
agreement between both codes, along with the low numerical 
diffusivity of FARGO, have validated FARGO as a very use- 
ful tool for studying in the embedded protoplanet problem. 
Since it is much faster than NIRVANA, it has been used 
to calculate the very long-term behaviour of an accreting 
protoplanet (see run F6 below). 

5.3 RH2D 

To obtain results using a third method we have employed the 
code RH2D (Kley 1989). This has been used previously in 
studies of disc-protoplanet interaction (Kley 1998, 1999a). It 
is a two-dimensional radiation hydrodynamics code. For the 
simulations presented here the radiation module is switched 
off, and all parts are solved explicitly. The code is based on 
the second order Van Leer (1977) advection algorithm and 
uses a staggered grid, with logarithmic spacing (a constant 
enlargement ratio of neighbouring grid cells) in radius and 
which is uniform in azimuth. As in FARGO, a fourth order 
Runge Kutta scheme was used to integrate the protoplanet 
orbit. This code has already been found to give similar re- 
sults to NIRVANA for protoplanet problems (Kley 1998). 
Accordingly we here use RH2D to study the orbital evolu- 
tion of an embedded planet under slightly different condi- 
tions to those adopted in the case of the other two numerical 
methods. 



6 NUMERICAL CALCULATIONS 

The main results of the numerical calculations are presented 
in table ^. Of particular significance are the migration times, 
T m i g , which are listed in the fifth column of table |l] in units of 
10 4 Po, and the estimated final masses of the planets, mfi na i, 
which are listed in the sixth column in units of the original 
planet mass mi n u. 

The values of r m i g were obtained by measuring the rate 
of change of the planet's semi-major axis at the end of the 



simulation and extrapolating forward in time assuming that 
this rate remains constant. Thus T m i g gives as estimate of 
the time required for the planet to migrate all the way to 
the central star. The values of mfi na i were obtained by mea- 
suring the accretion rate at the end of the simulation and 
extrapolating forward for a time r m i 9 . 

Thus m final provides an estimate of the mass the planet 
will attain on migrating all the way to the central star. It 
should be noted that while these values provide reasonable 
estimates of the true migration times and accreted masses, 
the assumption of there being a constant migration or ac- 
cretion rate is not strictly correct since these will tend to 
decrease as the planet mass continues to grow and mi- 
grate inwards. The only difference between runs Ni and Fi 
(i £ [1, 5]) is the code used to perform the calculations, with 
the Ni runs being performed with NIRVANA and the Fi runs 
with FARGO. The evolution equations of the system and the 
initial conditions are identical. Their results agree reason- 
ably well, except for the low resolution non accreting case, 
but as we shall see below, the low resolution we used for these 
runs is probably too coarse to give trustworthy results. From 
both the mid-resolution (N3 & F3) and high resolution (N5 
& F5) runs, we can see that FARGO gives slightly higher ac- 
cretion rates (and consistently longer migration times, since 
the inertia of the planet increases and, less importantly, most 
of the accreted material comes from the outer disk with a 
larger specific angular momentum). This higher accretion 
rate arises because FARGO has a smaller numerical diffu- 
sivity along the direction of orbital motion (Masset 1999). 
Associated with this is a more strongly peaked density pro- 
file around the planet. 

For the purposes of illustrating the main results of our 
simulations, we will now describe the results for one indi- 
vidual case. Following this we will compare the results of 
simulations in which the disc either did or did not have a 
tidally induced gap in which the planet orbits at the start 
of the calculation. We will then go on to compare the mi- 
gration rate of a planet which accretes gas from the disc at 
a maximal rate to that of a non accreting protoplanet. We 
then look at the effect of changing the initial surface density 
profile and disc aspect ratio. Following this, we compare the 
results obtained with NIRVANA and FARGO , and study 
the effects of changing the numerical resolution. 

6.1 An Illustrative Case 

The evolution of the protoplanet embedded in the disc for 
calculation N3 is shown in fig. [j], where grey-scale plots show- 
ing the surface density variation in the disc are presented, 
and the position of the planet is indicated by the small white 
circle (which has a radius equal to the Roche radius of the 
planet). Initially at time t — the disc surface density was 
unperturbed, but as the calculation proceeds the tidal force 
due to the planet strongly perturbs the disc, leading to the 
formation of trailing spiral shock waves. In particular an 
m = 2 spiral wave pattern may be observed. The transfer 
of angular momentum between the disc and the protoplanet 
leads to the formation of an annular gap, or surface den- 
sity depression, in the vicinity of the planet's orbit, which is 
cleared after about 200 orbits for a Jupiter mass planet. 

As the disc-planet system evolves, the inner disc is lost 
from the system since viscous evolution causes it to drain 
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Run 


Resolution 


Accretion 


Initial 


Tmig 


f i TLQ, I 




N r X N,f, 


on ? 


Gap ? 


xl0 4 P 




Nl 


50 x 147 


Yes 


No 


1.77 


4.8 


N2 


^0 v 1 A7 


No 


No 


1 07 




N3 


80 x 235 


Yes 


No 


1.00 


3.2 


N4 


80 x 235 


No 


No 


1.20 


1.0 


N5 


130 x 384 


Yes 


No 


0.84 


2.7 


N6 


130 x 384 


No 


No 


1.0 


1.0 


N7 


80 x 235 


Yes 


Yes 


0.45 


3.25 


N8 


80 x 235 


No 


Yes 


0.40 


1.0 


Fl 


50 x 148 


Yes 


No 


1.99 


4.77 


F2 


50 x 148 


No 


No 


1.75 


1.0 


F3 


80 x 236 


Yes 


No 


1.08 


3.85 


F4 


80 x 236 


No 


No 


1.33 


1.0 


F5 


130 x 384 


Yes 


No 


0.85 


3.27 


F6 


70 x 180 


Yes 


Yes 


2.03 


4.87 


Rl 


128 x 128 


(No) 


Gap 


1.45 


1.0 


R2 


128 x 128 


Yes 


Gap 


1.48 


4.63 


R3 


128 x 288 


Yes 


Gap 


1.36 


4.10 



Torques v's Time 



Table 1. The first column gives the run label (Ni NIRVANA, 
Fi FARGO, Ri RH2D), the second column gives the number of 
grid cells used, and the third column indicates whether the planet 
accretes from the disc. The fourth column indicates whether the 
calculation started with an initial gap in the disc. The fifth column 
gives the estimated migration time (in units of 10 4 Prj) and the 
sixth column provides an estimate of the final mass of the planet. 
Note that the estimated times are given at different evolutionary 
times of the models. 



through the open inner boundary. The tidal interaction be- 
tween the disc and protoplanet leads to the disc interior to 
r v exerting positive torques on the protoplanet, and the disc 
exterior to r v exerting negative torques on it. The loss of the 
inner disc leads to a reduction of the positive torque, so that 
the torque due to the outer disc becomes dominant. Conse- 
quently the planet undergoes inward orbital migration as its 
angular momentum is removed by the outer disc. In the final 
panel of fig. |l|, which corresponds to a time of t = 527 P , 
it may be observed that the semi-major axis of the planet's 
orbit is ~ 80 percent of it's original value, and that most of 
the inner disc has disappeared. As the calculation proceeds 
beyond this point the planet continues to spiral in towards 
the central star. 

The evolution of the torques acting on the protoplanet 
during the first 300 orbits are shown in fig. [], and indicate 
what the dominant contribution to the orbital migration is. 
Here, the torque per unit mass acting on the protoplanet due 
to the disc exterior to the orbital radius of the protoplanet 
r v (dotted line), interior to, or at, r p (dot-dashed line), and 
the indirect term in the gravitational force (dashed line) are 
plotted against time. Also plotted is the effective torque per 
unit mass that arises from the accretion of material whose 
specific angular momentum differs from that of the orbit 
(solid line). The first thing to note is that the accretion of 
gas from the disc has a negligible effect on the orbital evo- 
lution. Once the gap has formed the effective torque arising 
from accretion is only ~ 3 percent of the torque due to the 
outer disc. Second, we note that the torque arising from the 
outer disc material is consistently larger than that due to 




(Orbits) 



Figure 2. This figure shows the torques (per unit mass) acting 
on the protoplanet during the early stages of its evolution, in 
run N3. The torque due to the disc exterior to the planet orbital 
radius r p is shown by the dotted line, that due to the disc interior 
to or at r p is shown by the dot-dashed line, and the torque due 
to the disc's contribution to the indirect term of the gravitational 
potential is shown by the dashed line. The effective torque due to 
the accretion of disc gas is shown by the solid line. 



the inner disc material. The loss of the inner disc through 
the open inner boundary causes this disparity in the torques 
to grow, and ensures that the migration is always directed 
inwards. The contribution that the disc gravity makes to 
the indirect term is found to have a negligible effect on the 
orbital evolution. 

The evolution of the orbital radius of the protoplanet is 
shown by the solid line in the first panel of fig. ^, and that 
of the planet mass is shown by the solid line in the second 
panel. The evolution of the protoplanet undergoes a more 
rapid phase early on as the planet loses angular momentum 
to the disc material during the gap opening phase, leading to 
a faster period of orbital migration. This phase of evolution 
lasts for a few hundred orbits until the gap is opened, and 
also results in a very large mass accretion rate that almost 
doubles the mass of the protoplanet within t — 300 Po. 
Once the gap has been cleared, the evolution slows down 
and the planet spirals in towards the star at an almost con- 
stant rate while accreting gas from the disc at an almost 
steady accretion rate. The migration rate is observed to slow 
down slightly as the calculation proceeds since the increase 
in mass of the planet increases its inertia. The calculation 
was initiated with Mdo ~ 2m p , but as the planet migrates 
and accretes from the disc we arrive at a situation where 
m p > Mdo (see discussion in section [^). The orbital evo- 
lution of a massive planet orbiting in a disc that contains 
a smaller mass than the planet within a characteristic ra- 
dius r v has been studied by Ivanov, Papaloizou, & Polnarev 
(1999). In this physical regime the evolution is controlled 
by the viscous evolution of the disc and the inertia of the 
planet, with a more massive planet migrating more slowly, as 
described in section 2.2. We can test whether Eq. (H) agrees 



with our numerical calculations by examining the changes in 
planet mass, orbital radius, migration rate, and disc surface 
density at two different times, ti and £2, during the calcula- 
tion. If we allow for the reduction of the surface density in 
the outer disc due to accretion onto the planet, and by mass 



© 0000 RAS, MNRAS 000, 000-000 



10 Nelson, Papaloizou, Masset, Kley 



Orbital Radius v's Time Planet Mass v's Time 

1 .2 I 3.0 | 




second panel shows the increase in mass as the protoplanets accrete gas from the disc for the same calculations, with the same line styles 
as described above. 



flow through the gap and through the inner boundary, then 
we obtain the following expression relating the migration 
rates r p (t) 

r P (ti) = ( m p (t 2 )S(ti) \ 4/5 / r p (ti) \ 3/5 
f p (t 2 ) \m p (ti)T,(t2) J \r p (t 2 ) J 

for a disc with a constant value of v. Comparing the mi- 
gration rates in the numerical calculation at t\ — 1000 Po 
and ti = 2000 Po, we find reasonable agreement between 
the prediction of Eq. ( |22[ ) and the numerical results. The 
measured ratio of the migration rates gives r p (ti)/r p {t2) — 
1.5, whereas that predicted by Eq. (^2|) gives a value 
fp(ti)/r p (t2) — 1.4 with E taken as the azimuthally aver- 
aged value at the location of the 2:1 outer Lindblad reso- 
nance. 

The accretion rate decreases as the protoplanet mass 
increases because the disc is tidally truncated more effec- 
tively by a more massive planet (e.g. BCLNP) After a time 
of t ~ 2400 Po, the protoplanet has migrated to a radius 
of r ~ 0.6 and has accreted ~ 1.3 times its original mass. 
By extrapolating the migration and accretion rates forward 
in time, it is estimated that the planet will spiral into the 
central star after t ~ 1.0 x 10 4 Po, by which time it will 
have reached a mass mfi na i ~ 3.2 times its original mass. 

We note that the planet remains in an almost circular 
orbit throughout its evolution, and shows no sign of eccen- 
tricity growth, since the gap region still contains sufficient 
corotating material that is able to damp the eccentricity 
growth caused by the outer Lindblad resonances (Goldreich 
& Tremaine 1980, Artymowicz 1993a, b). 

The viscous evolution time of the disc is given by 
Eq. (^j) . At a radius of r = 1, this corresponds to an evolution 
time of ~ 10 4 Po, very similar to what is observed for the 
migration time of the protoplanets in all of our simulations. 

This confirms the idea described in section E.2 and in 



Lin & Papaloizou (1986) that giant protoplanets undergo- 
ing tidal interaction with a protostellar disc should migrate 
on a time controlled by the viscous evolution time of the 
disc when the interaction is sufficiently non linear to open 
up a gap, and when the mass of the planet is less than or 
comparable to the disc mass with which it gravitationally 
interacts. 

6.2 The Effects of an Initial Gap 

As well as performing calculations in which the initial disc 
was unperturbed, we also performed calculations in which 
the initial disc contained a tidally induced gap around the 
vicinity of the planet. Here we will concentrate on the cal- 
culations labelled as N3 and N7 in table |l|. 

The calculation labelled as N3 is for an accreting proto- 
planet embedded in an initially unpertur bed disc, and was 
described in detail in the previous section |6.l| The calcula- 
tion labelled as N7 in table ^ is for an accreting protoplanet 
initially embedded in a disc which has a tidally induced gap 
at time t = 0. This initial condition was obtained by running 
a calculation with a non accreting planet on a fixed circular 
orbit for t ~ 300 Po, until a clear gap was formed and the 
surface density in the gap region became steady. A reflecting 
inner boundary condition was employed during this phase in 
order that the disc mass was conserved. The results of this 
calculation were then used as the initial conditions of calcu- 
lation N7, but now with an accreting protoplanet which was 
able to undergo orbital evolution, and with an open inner 
boundary condition. 

The evolution of the orbital radius of the planet in cal- 
culation N7 is shown by the dotted line in the first panel of 
fig. [| and the evolution of the mass of the planet is shown 
by the dotted line in the second panel. In the case of the cal- 
culation N3, in which the planet is initially embedded in an 
unperturbed disc and must clear material in order to form a 
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gap, the clearing of that material leads to a period of more 
rapid migration. At the same time, this larger migration rate 
is slowed by the rapid accretion of gas, when the planet is 
deeply embedded, that leads almost to a doubling of the 
planet mass and inertia within a few hundred orbits. In the 
case of calculation N7, the planet initially resides within a 
gap, and so does not have to clear much material away from 
its vicinity during the early stages of its orbital evolution. 
It does, however, have an early period of more rapid accre- 
tion since it absorbs the material that is initially within its 
Roche lobe that accumulated there during the formation of 
the gap. This initially large accretion rate is augmented by 
the fact that an accreting planet helps to reduce the surface 
density of material in the gap region by accreting some of it. 
The similarity of the migration rates for the calculations N3 
and N7 during the first 1000 Po indicates that the effects of 
opening the gap for an initially fully embedded planet are 
almost entirely counter balanced by the planet mass growth 
and accretion of material with a higher specific angular mo- 
mentum, with this latter effect being almost negligible after 
~ 100 P (as shown in fig. ||). 

Looking at the the later stages of the evolution in the 
first panel of fig. ^, we notice that although the planet in 
calculation N3 migrates slightly faster to begin with, the 
planet in calculation N7 eventually migrates at a higher rate 
since its mass and inertia are smaller, with the orbital radii 
crossing over at t ~ 1500 Po- The expected migration time 
scale in calculation N7 is r m i g ~ 0.45 x 10 4 Po, with the 
final mass estimate being mfi na i ~ 3.2, indicating that the 
presence or not of an initial gap has only a relatively small 
effect on the final results. 



6.3 Comparison of the Evolution of an Accreting 
and Non Accreting Protoplanet 

Calculations were performed for both accreting and non ac- 
creting planets. In this section we will compare the results 
of calculations N3 and N4 in order to ascertain the effects 
of accretion on the evolution of the protoplanet. These cal- 
culations were both initiated with unperturbed disc models. 
For the accreting protoplanets, the accretion rate adopted is 
such that the e-folding time for the protoplanet to accrete 
mass within a distance of half of its Roche radius is given 
by Tacc = 3Po/(2-7r). This is a dynamical time scale so the 
situation corresponds essentially to a maximally accreting 
protoplanet (Kley 1999a). 

The case of a non accreting planet corresponds to the 
situation where the protoplanet has a lobe filling gaseous 
envelope in hydrostatic equilibrium with Kelvin-Helmholtz 
time longer than the migration time. Then the planet can 
accrete little mass during the migration stage. Such a situa- 
tion may occur if the envelope is built up while maintaining 
a low luminosity emitted from the central parts of the pro- 
toplanet. 

Calculation N3 was discussed in some detail in sec- 
Fig. ^ shows the evolution of the orbital separa- 



Orbital Radius v's Time 



tion 6.1 



tions for calculations N3 (solid line) and N4 (dotted line). 
It is apparent that the non accreting planet (N4) undergoes 
a significantly more rapid phase of migration as the initial 
gap is cleared during the first few hundred orbits, since the 
planet has to transfer a substantial amount of angular mo- 
mentum to the disc gas in order to clear the gap. The rapid 



1000 1500 
Time (Orbits) 



Figure 4. This figure shows the evolution of the orbital radius 
for the accreting protoplanet in run N3 (solid line) and for the 
non accreting protoplanet in run N4 (dotted line). 



accretion of mass (and some angular momentum) by the ac- 
creting planet (N3) helps to counteract the initial inward 
torques that arise when the gap is being cleared of material, 
as described in section |5,2| . 

As the calculations proceed, the migration rates in both 
cases slow down when the gap has been cleared of material. 
The planet then migrates inwards on approximately the vis- 
cous evolution time scale of the outer disc. Some additional 
slowing down of the migration rates occurs because the pro- 
toplanets interact with a smaller amount of disc mass as 
they migrate inwards. 

As the evolution time approaches t = 2000 Po, it is ap- 
parent that the migration rate of the non accreting planet 
is actually slower than that of the accreting planet, even 
though the accretion of material has increased the inertia 
of the accreting planet. The reason for this unexpected be- 
haviour is that the non accreting planet undergoes Roche 
lobe overflow. As material accumulates onto the non accret- 
ing planet, the Roche lobe eventually becomes filled with 
a hydrostatic atmosphere and no additional material may 
enter it. The continued flow of material from the outer disc 
onto the protoplanet then leads to circulation around it and 
Roche lobe overflow, such that material flows towards the 
central star. 

This material contains too much angular momentum to 
flow through the inner boundary directly, and instead fuels 
the inner disc. This inner disc then exerts a positive torque 
on the planet and reduces the rate at which it is able to 
migrate towards the central star. This process provides an 
efficient method of allowing material to flow across tidally 
induced gaps in accretion discs, and thus for the outer disc 
to feed material in to the inner disc which can continue to 
accrete onto the central star. This process will be the subject 
of a more detailed future study. 



6.4 A Long-term Evolution Run 

The run F6 is a low/mid-resolution run aimed at computing 
the long-term behaviour of the accretion/migration process. 
As indicated in table 0, N r = 70 and N4, = 180. Contrary 
to the other FARGO and NIRVANA runs, the inner bound- 
ary is located at r = 0.25 instead of r = 0.4, and the radial 
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Figure 5. This diagram shows the distance between the proto- 
planct and the primary as a function of time in run F6. At t = 
the gap is already cleared (it was cleared during 400 fixed circular 
orbits). 



Figure 6. This figure shows the temporal behaviour of the pro- 
toplanet mass as a function of time in run F6. We notice that the 
curve slope (i.e. the accretion rate) is not very peaked at t = 
since the gap was already present. The mass unit is one Jupiter 
mass. 



grid spacing is in geometric progression rather than uniform. 
Since we want to deal with the long-term behaviour of the 
protoplanet, we allow it to get close to the primary which 
is why we take a smaller inner boundary radius. Further- 
more, we want to track the accretion rate as accurately as 
possible. Taking the non uniform radial grid spacing ensures 
that the cells all have the same shape and that the accre- 
tion algorithm will not be biased accordingly. For this run we 
first clear the gap with the inner boundary open by evolving 
the system with the protoplanet orbit circular and fixed for 
400 periods. Then we start the accretion/migration process. 
Time t = is taken then. 

In this run, as distinct from the others, the frame is 
centered on the centre of mass of the system composed of the 
primary and the protoplanet. This is not an inertial frame, 
since that would need to be centered on the centre of mass 
of the primary, protoplanet and disc. But the indirect term 
arising from the acceleration of this frame is much smaller 
than the indirect term in the case of a frame centered on 
the primary. Furthermore, the material in the outer disk 
tends to orbit around the centre of mass of the primary and 
protoplanet (since the disk itself is not self-gravitating), so 
one can work with a rigid outer boundary and impose there a 
fixed Keplerian velocity, with no radial motion, which avoids 
inflow and outflow at that boundary. On the other hand, the 
material in the inner disk tends to orbit around the primary, 
so there is a mismatch there between the grid boundary 
and the gas orbits, which leads to a "vacuum-cleaner" effect 
which drains the inner disk faster than a frame centered on 
the primary would do. However, tests have shown that the 
flux of mass at the inner boundary is at most 10 to 20 % 
larger as a result of this effect. 

We present in figure |E] the evolution of the protoplanet- 
primary separation as a function of time. In figure ^ we show 
the total mass of the protoplanet as a function of time, and 
in table [j] we give the estimated migration time and final 
masses. These rates have been extrapolated from the time 
derivatives of the mass and orbital radius at time t = 7500. 
This corresponds roughly to the time at which the inner 
disk is lost. The fact that these results are in relatively good 
agreement with the previous ones, even though the extrapo- 



lation is performed after many more orbits in the case of run 
F6, shows that assuming that the migration and accretion 
rates reach constant values is a reasonable approximation, 
even though the curves show some residual deviation from 
linearity. 

We show in figure |^ a sequence of four surface density 
plots at times (t = 0, t = 300, t = 3000, t = 7000). As a 
consequence of the migration the gap radius decreases with 
time. We note as well the depletion of the inner disc (the 
viscous time scale at r = 0.7 is r v i sc = 7.8 ■ 10 3 orbits, and 
is even smaller at smaller radii). 

Fig. |^ shows the flow around the protoplanet at the end 
of the run. Even though the inner disk is strongly depleted, 
some Roche lobe overflow into it is indicated. The mean pro- 
file of the gap surface density at different times is displayed 
in fig. []. The gap deepens as the protoplanet mass increases 
with time. The depletion of the inner disc is also apparent. 

The mass lost through the inner boundary as a function 
of time is plotted in fig |l^. From the comparison of figs. ^ 
and [l^, one can see that the mass overflow flux and the 
mass accretion rate onto the planet have the same order of 
magnitude (the mass overflow flux is of the same order of 
magnitude as the mass outflow through the inner boundary 
since the mass lost at the inner boundary is bigger than the 
inner disk mass; in a totally stationary case these two rates 
would be strictly equal). 

The value of the ratio of the mass flow rate through the 
gap to the accretion rate ~ 1/3.5 obtained here is somewhat 
higher than that obtained from the simulation R2 where this 
ratio was about 1/7. However, the magnitude of the accre- 
tion rate is smaller here since the planet mass is larger. Also 
the surface density profiles and disk aspect ratios are differ- 
ent in these simulations. Given the differing physical situ- 
ation when the accretion rates were measured, we consider 
the agreement to be satisfactory. 
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Figure 7. Time sequence of the surface density for run F6. The outer boundary is at 2.5 and does not appear on the plots. The circle 
around the protoplanet in each case has a radius equal to the Roche radius. 



6.5 Runs with an Initial Surface Density Profile 
Using RH2D 

The code RH2D has been already been tested against NIR- 
VANA (Kley 1998) and found to give very similar results 
so we shall not give results of additional tests here. Instead 
we present runs Ri (i — 1,2,3) which use different initial 
conditions incorporating a surface density profile that is not 
constant, as outlined at the end of section 4.6, and a slightly 
higher value of the disc aspect ratio H/r — 0.05. Also, it 
should be noted that these calculations are initiated with a 
gap already in existence around the position of the proto- 
planet. 

The main results concerning the mass accretion by the 
protoplanet are displayed in Fig. The lower panel shows 



the evolution of the protoplanet mass and the mass lost 
through the inner boundary for the models R2 and R3 which 
have maximal accretion. These models have the same physi- 
cal setup but different numerical resolutions. Note that only 
the mass added to the planet during the evolution is dis- 
played. The total planet mass is obtained by adding lMj to 
the quoted values, i.e. at the end of the run R2 at t ~ 3300 
the planet has reached about 2M.j. This value is smaller than 
that obtained by the other codes (see for example run F6), 
which is a result of the surface density near the planet being 
a factor of about three smaller in these runs. Also shown in 
the lower panel is the mass lost through the inner boundary 
(dotted and dashed dotted lines) which is assumed to have 
been accreted by the star. The mass of the planet rises more 
slowly initially, because of the initial gap imposed. Only 
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Figure 8. Details of the flow in the neighbourhood of the pro- 
toplanet near the end of the run F6 (t = 7000Prj)- The velocity 
field is represented in the frame corotating with the planet. We 
clearly see the gap, the wakes of the protoplanet in the inner 
and outer disk, the extremities of the horseshoe orbits in the 
gap, the two X-points in the velocity field, which correspond to 
the Lagrange points Li and L2. The circle is centered on the 
protoplanet and has a radius equal to the Roche radius of the 
protoplanet. Because the actual potential felt by the disk mate- 
rial is smoothed, the X-points lie slightly inside this circle. One 
can see that the density peak around the planet is not elongated 
along the orbit, as usually observed in fixed frame normal advee- 
tion scheme runs. Because of the geometric grid spacing, all the 
cells have the same shape and are "as square as possible" , i.e. 
log(r out/Tin) ~ AT r log(l + 2-k/N^). 
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Figure 10. Mass lost at the inner boundary as a function of time 
in run F6. At t = this mass is not zero since the run was begun 
with a non-accreting protoplanet on a fixed circular orbit with an 
open inner boundary in order to generate an initial gap. 
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Figure 11. Mass evolution for models R2 and R3: The lower 
panel shows the evolution of the mass of the planet and the mass 
lost from the disc (accreted by the star). The top panel shows 
the corresponding accretion rates in units of Mj/Pq. Dark solid 
and dotted lines refer to model R2, while the lighter dashed and 
dashed-dotted lines refer to the higher resolution model R3. 



Figure 9. Disk surface density profile at different times for run F6 
(t = : solid line, t = 690 : dotted line, t = 1830 : dashed line and 
t = 3550 : dot-dashed line). These plots represent the azimuthal 
average of the surface density, and hence take the protoplanet 
wake into account. The residual surface density which would be 
obtained by omitting the wakes would be much smaller. 



after the full development of the quasi-stationary flow (at 
t w 300) does the mass accumulation rate onto the planet 
become larger than the mass loss rate through the inner 
boundary. This is demonstrated in the top panel where the 
mass accretion rates in units of Mj/Po onto the planet and 
star are shown. After the mass contained initially within 
r = 0.25 and 0.4, Mdo = 0.37Mj is consumed (primarily by 
the star), the mass accretion rate onto the star settles to the 
very small constant value Mi os t = 2.86 x 10 — Mj/Po which 



© 0000 RAS, MNRAS 000, 000-000 



Migration of protoplanets 15 
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Figure 12. Radial evolution for models Hi: The lower panel 
shows the evolution of the semi major axes of the planet for three 
different models. Solid black line R2, dotted grey line the non- 
accreting model Rl, and dashed light grey line the high resolution 
model R3. The top panel shows the corresponding derivative r p 
in units of ao/Po versus time. The diamonds represent the analyt- 
ical approximation according to Ivanov, Papaloizou, & Polnarev 
(1999) as given in Eq. M for model R2. 



proximation according to Ivanov, Papaloizou, & Polnarev 
(1999) as given in Eq. |^ with an arbitrary normalisation 
factor which can be adjusted to match the numerically ob- 
tained data. For simplicity here the maximum of the az- 
imuthally averaged surface density E(r) outside of the planet 
was taken for model R2 as the value for the surface density 
in (^|). Clearly the formula gives an excellent approximation 
to the actual evolution of the radius of the planet. 

6.6 Comparison between NIRVANA and FARGO 

In this section we present the results of calculations that 
were performed with both NIRVANA and FARGO in or- 
der to check that the results that we have presented are 
reproducible when using independent numerical codes. The 
calculations that we will compare are N3 performed using 
NIRVANA and F3 performed using FARGO. As described in 
table |l[ these calculations are for an accreting Jupiter mass 
planet initially embedded in an unperturbed disc model. 
For comparison purposes, the orbital evolution of the pro- 
toplanet is plotted in the first panel of fig. [l3|, with calcu- 
lation N3 being represented by the solid line and F3 be- 
ing represented by the dotted line. Similarly, the evolution 
of the protoplanet mass is presented in the second panel 
of this figure with the same line styles being used to rep- 
resent the calculations. It is apparent that the results are 
very closely matched in terms of the orbital evolution rate 
and reasonably well matched in terms of the mass accretion 
rates. We clearly see the trend, also seen in the results given 
in table [l], that the FARGO algorithm gives higher accre- 
tion rates and hence slightly slower migration. This may be 
related to the fact that the distribution of matter close to 
the protoplanet is subject to less numerical diffusion and so 
has less azimuthal elongation in FARGO. 



may be identified with the rate of mass overflow across the 
gap. Near the end of the computation at t = 3000 this rate 
is about 6.5 times lower than the mass accretion rate by the 
planet which is in rather good agreement with the results 
quoted in Kley (1999a) but somewhat lower than the rates 
obtained with FARGO as outlined above. 

In Fig. [l^ the evolution of the distance of the planet 
from the star is shown for all models R1-R3. Again, due to 
the smaller mass contained in the physical domain the or- 
bital decay rate is smaller than in the previous models (Ni 
and Fi). In the comparison model labeled 'constant mass' 
(Rl, dotted line) the mass of the planet remained constant 
but mass was nevertheless removed from the Roche-lobe in 



the manner described in section 4.5 (also see Kley, 1999a). 
This model has a smaller decay rate than the complete mod- 
els (R2 solid line, and R3, light dashed line) in which the 
accreted mass is added to the planet mass. All the models 
show initially a fluctuation in r p (t) which indicates an eccen- 
tricity growth to less than about 0.02 which damps out at 
later times. For the high resolution model (R3, light dashed 
line) the radial evolution follows closely the low resolution 
model, the migration rate seems to be marginally lower. 
This is consistent with the results obtained with NIRVANA 
and FARGO which indicate that lower resolution runs have 
larger accretion rates and slower decay rates (see Figs, [lo] 
and|l|). 

The diamonds in Fig. ^ represent the analytical ap- 



6.7 Numerical Resolution 

6. 7. 1 Non accreting planets 

In this section we present the results of calculations that 
were performed using different resolutions. We first concen- 
trate on comparing three calculations in which the orbital 
evolution of a non accreting planet was studied, namely cal- 
culations N2, N4, and N6. 

The evolution of the planet orbit radius for these cal- 
culations is presented in fig. [l4|, with calculation N2 being 
represented by the dotted line, calculation N4 by the solid 
line, and calculation N6 by the dashed line. The agreement 
between N2 and the other calculations appears to be the 
worst, which is to be expected since this is the lowest res- 
olution simulation that we performed and this is probably 
too low to fully resolve small scale structures in the vicin- 
ity of the planet. The two calculations N4 and N6 on the 
other hand show extremely good agreement in their orbital 
migration rates at later times, though there is a small offset 
in the orbital radius at any given time due to calculation N6 
experiencing more rapid migration during the gap clearing 
stage of the calculation. 

This increased rate of orbital migration for the higher 
resolution calculation during the gap clearing stage proba- 
bly arises because of its ability to resolve density waves with 
higher azimuthal mode numbers that are located close to the 
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F3 by the dotted line. The second panel shows the increase in mass as the planet accretes gas from the disc, with the same line styles 
representing the two calculations. 
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Figure 14. This figure shows the evolution of the orbital ra- 
dius for the calculations N2 (dotted line), N4 (solid line), and N6 
(dashed line). The resolution for each calculation is described in 
table |. 

planet when it is embedded. Once the gap has been cleared, 
however, these non axisymmetric structures do not provide 
a significant contribution to the tidal torque acting on the 
planet, so that the migration rates then become approxi- 
mately equal. The close agreement between the migration 
rates of calculations N4 and N6 indicate that our calcula- 
tions have essentially reached convergence in their results. 

6.7.2 Accreting planets 

The calculations for accreting planets initially embedded in 
unperturbed accretion discs that were performed with differ- 
ent resolutions are shown in figs. [L5| and jH]. The first panel 
in fig. [l5| shows the orbital radius as a function of time for 
the calculations Nl (dotted line), N3 (solid line), and N5 



(dashed line). The second panel shows the evolution of the 
planet mass, with the same line styles as described for the 
first panel. 

Fig. ^ gives the corresponding results for the runs Fl, 
F3 and F5. It is apparent that the lower resolution runs have 
a larger accretion rate. This leads to smaller migration rates 
because the protoplanets have larger masses. Nonetheless, 
the migration times obtained for these runs all indicate an 
orbital decay time of r m i g ~ 10 4 Po, as indicated in table |l[ 
This is in agreement with the idea that the orbits of giant 
planets will evolve on the viscous diffusion time of the pro- 
tostellar disc. The results of NIRVANA and FARGO are in 
good agreement with respect to migration rates and final 
planetary masses. The predicted final planetary masses in 
the calculations N3 and N5 are 3.2Mj and 2.7Mj while for 
the calculations F3 and F5 they are 3.85M./ and 3.25Mj 
respectively. These agree well with the estimated masses of 
a number of the recently discovered closely orbiting extra- 
solar giant planets (Marcy, Cochran, & Mayor 1999, Marcy 
& Butler 1998). 



7 DISCUSSION AND CONCLUSION 

In this paper we have studied the interaction of a proto- 
planet of lMj initially with a gaseous disc whose aspect 
ratio and kinematic viscosity are those expected for a mini- 
mum mass solar nebula. This had characteristically 2Mj in- 
terior to the initial circular orbit radius of the protoplanet. 
The problem was studied with three independent hydrody- 
namic codes, NIRVANA, FARGO and RH2D. These were 
found to give consistent results when compared. FARGO 
had the additional advantage that, on account of the fast 
advection scheme employed, the evolution could be followed 
for a much longer time. 

A general result of the simulations was that the direc- 
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Figure 15. The first panel shows the evolution of the orbital radius for the calculations Nl (dotted line), N3 (solid line), and N5 (dashed 
line). The evolution of the planet mass is shown in the second panel, with the line styles being the same as above. The resolution for 
each calculation is described in table ll. 
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Figure 16. Same as figure [Ti] for FARGO results Fl (dotted line), F3 (solid line) and F5 (dashed line). 



tion of the orbital migration was always inwards and such 
that the protoplanet reached the central star in a near circu- 
lar orbit after a time of ~ 10 initial orbital periods, which 
is characteristically the viscous time scale at the initial or- 
bital radius. This was found to be independent of whether 
the protoplanet was allowed to accrete mass or not, or the 
surface density profile in the disc. The tendency to migrate 
inwards was assisted by the disappearance of the inner disc 
through accretion onto the central star. When the proto- 
planet was allowed to accrete at a near maximal rate the 



mass was found to increase to about 4Mj as it reached the 
central star. Because of deep gap formation and lower ac- 
cretion rates for the larger masses (also see BCLNP) it is 
difficult to exceed this mass in the kind of simulations pre- 
sented here. An additional calculation has been performed 
with an initial planet mass of 3Mj, but whose results are 
not described here in detail. This calculation resulted in a 
similar migration time of ~ 10 4 Po, and an estimated final 
mass for the protoplanet of 4.8 Mj, indicating the difficulty 
of forming giant planets with masses greater than about 5 



© 0000 RAS, MNRAS 000, 000-000 



18 Nelson, Papaloizou, Masset, Kley 



Mj before they have migrated to the centre. It would ap- 
pear that the masses estimated for a number of close orbit- 
ing giant planets (Marcy, Cochran, & Mayor 1999, Marcy & 
Butler 1998) as well as their inward orbital migration can be 
explained by the operation of the processes considered here 
during the late stages of giant planet formation. 

Several important issues, however, remain to be re- 
solved. The inward migration time is shorter than previously 
estimated time scales ~ 10 6 yr for the formation of a Jovian 
mass protoplanet (e.g. Bodenheimer & Pollack 1986; Pollack 
et al. 1996; Papaloizou & Terquem 1999). This is suggestive 
that either the planet formation process may be speeded up 
by the earlier merger of cores undergoing type I migration, 
and/or there may be regions in the disc where the viscos- 
ity is very small, so halting type II migration, perhaps due 
to inadequate ionization for MHD instabilities to operate. 
Additional planets embedded in the disc alter the density 
structure and consequently the torque balance which may 
result in a halting of the migration process (Kley 1999b). 
Many body processes such as gravitational scattering of pro- 
toplanets may also operate to move them to different orbital 
locations in the disc (e.g. Weidenschilling & Marzari 1996). 

An additional possibility is that in some cases giant 
planet formation occurs at substantially larger distances 
from the host star than have hitherto been given serious 
consideration. For example, a planet forming at a radius 
~ 20 AU will have a migration time of > 10 yr, which is 
now within the estimated range of lifetimes of protostellar 
discs around T Tauri stars (i.e. 10 6 - 10 7 yr). It is possible 
that inward migration of these protoplanets may simply be 
halted by the eventual dissipation of the disc at the end of 
the T Tauri stage. The final orbital positions of these plan- 
ets will then be determined by the initial radius at which 
the planets were formed and the age at which the T Tauri 
phase ends. In this scenario, planets that initially start to 
form closer in towards the central star (e.g. at 5 AU) will mi- 
grate inwards and will become 'hot Jupiters', whereas those 
planets that form further away stand a much greater chance 
of being at intermediate distances from their host stars when 
orbital migration is halted by the disappearance of the disc. 

Another issue is that type II migration in a viscous disc 
as considered here tends to cause Jovian mass protoplanets 
to merge with their central star on a time scale short com- 
pared to the lifetime of protostellar discs. Thus a process for 
halting the migration is required. This may occur through 
the termination of the inner disc due to a magnctospheric 
cavity (Lin, Bodenheimer, & Richardson 1996). 

The calculations presented here make a number of im- 
portant assumptions that may have some bearing on the fi- 
nal results obtained. By using a locally isothermal equation 
of state, we tacitly assume that any heat generated by the 
spiral shock waves is immediately radiated from the system. 
This may not be an accurate description of the thermody- 
namics, and if radiative processes operate on a time scale 
that is longer than the orbital time scale then some thicken- 
ing of the disc may result. In addition, we assume that the 
turbulent viscosity can be simply modeled using the Navier- 
Stokes equation, when in reality it should arise naturally 
from MHD instabilities. These, and other assumptions can 
only be addressed by performing global simulations that in- 
clude radiative and MHD processes in three dimensions. We 



hope in the near future to be able to address some of these 
outstanding issues. 
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